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When  using  a  numerical  method  of  approximating  the  value  of 
b 

the  definite  integral  5  f( x)dx,  it  ia  of  course  very  important 

a 

that  one  knows  the  character  of  t.he  integrand)  does  f(x)  or  its 
derivative  have  any  singularities  in  the  closed  interval  [a,b], 
and,  if  so,  what  type  of  singularities  are  present?  The  answers 
to  these  questions  will  determine  the  type  of  quadrature  formula 
which  should  be  used. 

In  the  case  of  a  one-dimensional  integral  with  a  reasonably 
simple  integrand,  one  can  usually  determine  fairly  easily  whether 
or  not  a  singularity  is  present  and  -  perhaps  with  more  difficulty 
the  character  of  the  singularity.  Sometimes  the  singularity  may 
be  removed  by  a  change  of  variables,  so  that  a  standard  non¬ 
singular  type  quadrature  formula  may  be  used.  In  other  cases  one 
uses  a  generalized  Gauss  formula  (or  similar  formula  with  other 
spacing)  which  has  been  derived  using  a  proper  weighting  function 
for  the  singularity  in  question. 

In  the  case  of  two-dimensional  integrals,  it  is  easier  to  be 
misled  by  one's  intuition.  With  the  advent  of  high-speed  digital 
computers,  it  has  become  the  tendency  to  ask  for  general  computer 
programs  which  will  integrate  any  "reasonable"  function.  Multiple 
integrals  are  usually  treated  as  repeated  simple  integrals,  (so 
that  one-dimensional  quadrature  formulas  are  used  two  or  more 
times).  It  is  the  purpose  of  these  notes  to  give  some  simple 
examples  of  multiple  integrals  which  show  some  of  the  difficulties 
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which  may  arise.  Also  some  formulas  will  be  given  to  use  on  singular 
integrals  of  the  types 

5h  jh  d)dY 

0  0  777? 

and 

h  h  r? — ? 

3  5  In  Jx  +  y  F(x,y)dxdy  f 

0  0 

which  occur  from  time  to  time  in  physical  problems.  In  the  three 
appendices,  the  formulas  discussed  in  the  text  are  collected  for 
ready  reference.  No  formal  expressions  for  error  terms  associated 
with  the  formulas  of  Appendix  II  an!  Appendix  III  are  given. 

Si.  Examples  in  which  the  integrand  is  non-singular. 

Example  1:  Consider  5  3  F(x,y)dxdy,  where  F(x,y)  =  1  and 

A  2 

where  A  is  the  first  quadrant  area  enclosed  by  the  curve  y  =  1  -  x  . 

(See  Figure  l) . 


Figure  1 
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The  multiple  integral  may  be  written  as  an  iterated  integral  in 
either  of  two  ways: 


(1.1) 


(1.2) 


With  either  form,  one  can  easily  perform  the  integrations  analytically 
and  arrive  at  the  exact  answer  I  =  2/3.  Since  F(x,y)  is  such  a  smooth 
function  (identically  l) ,  one  might  suspect  at  first  sight  that  (l.l)  or 
(1.2)  could  be  evaluated  quite  well  by  use  of  Simpson's  rule.  This  was 
actually  programmed  on  a  machine.  In  both  cases  the  inner  integrals 
were  evaluated  by  using  a  5 -point  Simpson  formula.  The  second  integra¬ 
tion  was  also  performed  using  Simpson's  rule,  and  it  was  tried  several 
times  using  different  integration  step  lengths.  The  results  are  given 
in  Table  1. 


Number  of  Points  Used 

In  Outer  Integration 

Numerical  Result 

Numerical  Result 

Formula 

Problem  (1.1) 

Problem  (1.2) 

5 

.65652626 

. 66666668 

9 

.66307927 

.66666667 

17 

.66539809 

.66666676 

33 

.66621795 

.66666672 

65 

.66650784 

. 66666688 

129 

.66660944 

.66666813 

257 

. 66664446 

.66666923 

Table  I:  Results  (using  Simpson's  Rule)  to  evaluate  integral  in  (l.l) 
and  (1.2) 


Comments  on  the  results ;  It  is  seen  that  for  problem  (1.2)  Simpson's 
rule  gives  very  close  to  the  correct  answer  2/3-  As  more  points  are 
taken,  the  results  do  get  poorer,  but  this  is  caused  entirely  by  round- 
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off  error. 

In  the  case  of  problem  (1.1),  Simpson's  rule  gives  results  which 
are  quite  poor.  The  results  do  get  better  as  more  points  are  taken, 
but  if  five-place  accuracy  were  desired,  many  more  than  the  257  *  5  =  1275 
points  would  be  required;  it  is  obvious  by  looking  at  the  differences, 
caused  by  round-off  in  problem  (1.2)  that  round-off  errors  would  also 
"take  over"  in  problem  (l.l)  before  five-figure  accuracy  could  be 
attained . 


With  such  a  simple  example,  it  is  quite  readily  seen  why  the  order 
of  integration  makes  so  much  difference.  If  the  inner  integration 
were  done  analytically,  equation  (l.l)  would  become 


1  _ 

I  =  J  -7l  -  y  dy 
0 

while  (1.2)  would  become 
1  2 

1  =  5  (1  -  x*)dx. 
0 


(1.3) 


(1.4) 


Theoretically,  (i.e.,  except  for  round-off),  Simpson's  rule  should 

give  exact  results  for  the  problem  of  integrating  the  polynomial 
2 

G(x)  =  1  -  x  over  the  interval  [0,1].  On  the  other  ha. i ,  the 
function  H(y)  =  «/l  -  y  may  not  be  approximated  very  accurately  by 
portions  of  cubic  polynomials  (because  it  has  an  infinite  slope  at 
y  =  1) .  Simpson's  rule  should  not  be  used  for  such  a  quadrature; 
instead  a  special  formula  should  be  used  which  takes  into  account 
the  nature  of  the  singularity  in  H'(y). 
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Example  2.  Consider  the  Integral 


j  =  s'  i^J7^=*y?T7j 
-1  -JZZ?IJ?77  ' 


dtdx. 


(1.5) 


Geometrically,  the  symmetric,  non-negative  integrand  is  being  integrated 

over  a  circle  of  radius  1.  The  integration  is  easily  performed 

analytically,  and  it  is  found  that  J  =  tt.  If  one  considers  this  as 

a  repeated  simple  integral,  and  if  a  twelve-point  Gauss  quadrature 

formula  is  used  in  the  intervals  [-  Jl  -  x?,  Jl  -  xf],J  =  1,2,..., 12 

J  J 

for  the  inner  integrations  and  in  the  interval  [-  1,1]  for  the  outer 
integration  (so  that  the  integrand  has  been  evaluated  at  144  points  in 
all),  the  result  3-1440  is  obtained;  the  error  is  +  .0024.  If  thirteen- 
point  Gauss  formulas  are  used  (so  the  integrand  has  been  evaluated  169 
times  in  all),  the  result  3.1379  is  obtained;  the  error  is  -  .0037. 

One  might  surmise  that  a  reason  for  the  very  inaccurate  results 
might  be  that  the  integration  with  respect  to  t  has  not  been  done 
very  accurately  when  x  is  close  to  zero.  For  example,  when  x  =  0, 
it  is  desired  to  integrate  | 2 t |  which  has  a  discontinuity  in  the 
derivative  when  t  ~  0.  Thus,  one  expects  better  results  if  he  evaluates 
numerically 
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If  six-point  Gauss  formulas  are  used  for  the  integrals  (1.6) , 
(equivalent  again  to  the  use  of  144  points  over  the  whole  circle) ,  the 
final  result  is  3.1437.  The  error  is  now  +.0021  (as  opposed  to  +.0024 
for  the  144  point  case  before) ,  so  the  error  is  still  more  than  expected 
for  the  integration  of  ’’smooth"  functions. 

Much  more  accurate  results  are  obtained  if  one  uses  the  ordinary 
Gauss  quadrature  formulas  for  performing  the  t  integration  and  then 
uses  the  generalized  Gauss  formula 

S  Jl  -  x2  f(x)dx  =  Hjf ( x^)  +  H2f(x2)  +•••+  H12f(x12)  (1.7) 

for  the  outer  quadrature.  Here  x^  =  cos  and  sin^ 

for  i  =  1,2,3, ... ,12.  The  final  result  using  this  formula  is 
J  -  3.141630.  The  error  is  only  +.000037. 

Comments  on  the  results:  It  is  often  stated  that  if  a  double  integration 
problem  is  symmetric  with  respect  to  the  two  independent  variables,  one 
should  "treat  both  variables  the  same".  Then  why  does  the  use  of  ordinary 
Gauss  integration  in  both  directions  yield  such  poor  results,  and  why 
are  results  improved  so  much  by  using. a  different  formula  for  the  x 
Integration?  The  answer  presents  itself  when  the  t  integration  is 
actually  performed  analytically.  The  result  is 
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1  / - 2 

Thus  we  now  wish  to  evaluate  25  vl  -  x  dx,  the  integrand  of  which 

-1 

has  infinite  slopes  at  x  =  +  1.  The  twelve-point  ordinary  Gauss 

formula  gave  better  than  four-place  accuracy  when  the  t  integrations 

were  performed,  and  the  results  were  of  course  good  approximations  to 

JT77\  The  generalized  Gauss  formula  will  evaluate  5  \/l  -  x2  dx 

-1 

exactly,  so  it  does  a  good  job  of  integrating  the  "good  approximation" 

to  Jl  -  x2. 

Hence,  for  a  symmetric  problem,  it  is  often  not  really  possible 
or  desirable  to  treat  both  variables  "the  same"  when  multiple  inte¬ 
gration  is  performed  by  iterated  use  of  one-dimensional  formulas. 

The  choice  of  a  particular' variable,  with  respect  to  which  one 
integrates  first,  ruins  the  symmetry  of  the  problem.  (These  remarks 
do  not  apply  if  the  region  of  integration  is  a  rectangle.) 

a  a 


52.  Singular  integrals  of  the  type  1=5  5 


—a  —a  /  2  2 

Jx  +  y 


dxdy . 


(2.1) 


Here  it  will  be  assumed  that  F(x,y)  is  analytic  in  the  square 
-  a  <  x,  y  <  a  and  that  F(0,0)  /  0.  The  difficulties  illustrated  in 
§1  with  problems  in  which  the  original  integrand  was  always  finite  and 
continuous  should  lead  one  to  expect  that  even  more  care  is  necessary 
in  the  use  of  numerical  methods  when  the  integrand  is  infinite  at 
some  point. 


Theoretically,  if  y  /  0,  the  integration  with  respect  to  x 
could  be  carried  out  with  a  standard  non-singular  formula  like  Simpson's 
rule.  This  would  give  an  approximation  for 
a 


G(y)  =  5 


dx. 


x  +  y 
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Since  F(x,y)  is  analytic  in  -  a  ^  x  <  a,  -  a  <  y  <  a,  it  may  be 
expanded  in  a  Taylor  series  (with  respect  to  x) ,  so  that 

2  2 

F(x,y)  -  F(0,y)  -  xFx(0,y)  =  fj-  0^(0, y)  +  |  F^O^y)  +•••}  =  §j-  H(x,y)  , 


where  H(x,y)  will  also  be  analytic  in  the  square.  Assuming  s'fcill 


that  y  /  0,  G(y)  may  then  be  written 


G(y)  =  F(0,y)  J 


+  F  ( 0,y)  5 


-a  /  2  ^  2 

n/x  +  y 


xdx  +  1_  ja  x2H(x,y)  d 

-a  12,  2  21  -a  /  2  ,  2 

Vx  +  y  vx  +  y 


(2.2) 


G(y)  =  F(0,y)  In 


■s/y2  +  a2  +  a  +  _1 
21 


i  5a 

-  JJT7 


(2.3) 


(since  the  second  integral  in  (2.2)  is  zero  because  the  integrand  is 
an  odd  function  of  x)  . 


The  integral  in  (2.3)  is  obviously  a  continuous  function  of  y 
which  will  approach  a  finite  limit  when  y->0.  The  first  term  in 
(2.3)  is  also  a  continuous  function  of  y  if  y  /  0.  But  at  y  =  0 
this  term  is  logarithmically  singular.  That  is,  (since  F(0,0)  =  constant/  0)  , 


=  -  2F(0,0) . 
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Thus  G(y)  may  be  written  in  the  form  g(y)  ln|y|,  where  g(y)  will 
remain  finite  throughout  [-  a, a]  if  0  <  a  <  1.  The  problem  of  per¬ 
forming  the  cubature  of  (2.1)  has  now  been  reduced  to  that  of  performing 
the  two  quadratures 


sW  *  I5M  d*  (*•*> 

lyl  -»  J7 77 


and 

a 

1=5  ln|y|  g(y)dy.  (2.5) 

-a 


Theoretically,  for  example,  the  standard  Gauss  formula  or  Simpson's 
rule  could  be  used  to  approximate  (2.4)  when  y  f  0.  A  generalized 
Gauss  formula  with  weighting  function  In  |y |  may  then  be  used  to 
approximate  (2.'5).  (For  such  a  formula,  see  Appendix  i)  .  Methods  of 
this  type  have  been  used  successfully,  but  here  also  it  is  easy  to 
overlook  some  difficulty  which  will  cause  large  errors  in  the  result. 
Consider  the  special  example: 


Example  3: 


1 

5 


o 


4 


+  x2  +  y 

Jx2  +  y2 


2 

-dxdy . 


This  integral  is  essentially  of  the  type  (2.1).  Since  the  integrand 
is  symmetric  in  x  and  y  the  two  lower  limits  were  taken  as  zero. 
Equations  (2.5)  and  (2.4)  for  this  problem  would  be 


1 

Io  =  5  ln|y|  g(y)dy 
0 


(2.6) 
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where 


The  four-point  generalized  Gauss  formula 

I3  ~  H]g(yi)  +  H2g^y2^  +  H3g(y3)  +  H4g(y4') 

where  H  =  -  .383464068  y1  =  .041448480 

H2  =  -  .386875318  y2  =  .245274914 

H3  =  -  .190435127  y3  =  .556165454 

H4  =  -  .039225487  '  y4  =  .848982395, 

was  used  to  evaluate  the  integral  in  (2.6).  (The  necessary  values' of 
g(y)  were  obtained  previously  by  using  an  ordinary  four-point  Gauss 
formula  applied  "to  the  integral  of  (2.7))  ..  The  resulting  approximate 
value  for  I3  was  7.58863*  The  integral  I3  may  of  course  be 
evaluated  analytically,  and  the. result  is 

I  =  |{25  ln(  1  +J2)  +  J2}  =  7.816186. 

The  approximate  result  was  too  low  by  about  0.23*  If  we  had  not 
worried  abput  the  singularity  at  all,  and  had  merely  used'  the  regular 
four-point  Gauss  formula  for  both  the  x  and  y  integrations,  the 
result  would  have  been  7.7219,  which  is  too  low  by  only  about  0.09. 

Comments  on  the  results:  The  question,  of  course,  is  why  better 

results  are  obtained  when  one  ignores  the  presence  of  the  singularity 
completely  than  if  formulas  (2.6)  and  (2.7)  are  used.  The  main  reason 
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is  that  in”taking  care  ofn  the  singularity  at  y  =  0,  a  new  singu¬ 
larity  was  introduced  at  y  =  1.  When  using  a  generalized  Gauss  formula 
to  evaluate  (2.6),  it  is  assumed  that  g(y)  may  be  approximated 
reasonably  well  by  a  polynomial.  However,  equation  (2.7)  shows  that 
g(l)  =  °°,  so  the  generalized  Gauss  formula  should  not  have  been  used 
to  evaluate  (2.6) j  to  use  it  was  a  definite  error  in  reasoning.  The 
example  was  given  here  however  to  emphasize  that  such  errors  will 
often  not  show  up  when  the  problem  is  placed  on  the  machine  unless 
the  exact  answer  is  known. 


The  integral  could  Have  been  approximated  better  by  treating 

a  smaller  region  about  the  singularity  as  a  singular  problem  and  by 
using  standard  quadrature  formulas  in  the  region  away  from  the  singu¬ 
larity.  This'  correct  use  of  formulas  analogous  to  (2.6)  and  (2.7)  is 
shown  in  the  evaluation  of  the  integrand  of  Example  3  over  a  smaller 
region:  ' 


Example  4: 


±  ±  ^  2  2 
I ,  =  5  5  - - - dxdy . 

0  0  •  J7T7 


This  integral  may  be  evaluated  exactly  and  the  result  is 

1 4  =  i^v/2  +  385  ln(  1  +  J2)  ]  §  1.7747034- 

If  the  singularity  is  ignored  and  an  ‘ordinary  four-point  Gauss 
formula  is  used ,• for ■ both  x  and  y  integrations,  the  resulting  value 
is  ■  1.7511  which  is  too  low.  by  about  .0236. 
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If  one  desires  to  perform  the  cubature  as  an  iterated  quadrature 
(but  making  use  of  the  fact  that  the  integrand  determined  by  the  inner 
integral  has  a  logarithmic  singularity  at  y  =  0) ,  the  specific  formula 
will  be 

1/  =  S  In  y  g(y)dy,  (2.8) 

4  A 


where 


5 

o 


.  .  2  ,  2 
4  +  x  +  y.  dx. 
nr— 2 

Jx  +  y 


(2.9) 


The  generalized  four-point  Gauss  formula  needed  for  this  problem  is  of 
the  form 


5 

o 


In  y  g(y)dy  =  H1g(y1)  +  H2g(y2)  +  H^gty^)  +  »4g(y4) . 


The  abscissas  y^  and  the  weights  are  given  in  the  table  of 

Appendix  I.  Using  this  formula  (and  using  ordinary  four-point  Gauss 
quadrature  formulas  for  approximating  each  functional  value  g(y^))> 
one  obtains  the  value  1.7728  for  I which  is  only  in  error  by 
about  0.0019. 


Comments  on  the  results  t  The  fact  that  the  error  0.0019  (incurred 
when  the  presence  of  the  singularity  is  considered)  is  only  about 
eight  percent  of  that  incurred  when  the  singularity  is  ignored  is 
encouraging.  However,  the  more  accurate  method  is  also  subject  to 
criticism.  Should  the  ordinary  Gauss  formula  be  used  to  evaluate 
(2.9)  for  the  various  desired  values  of  y?  It  is  true  that  for 
y  /  0  there  is  no  singular  point  in  0  <  x  <  but  as  y  becomes 
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very  small  the  function 

.,2.2 

4  +  x  +  y 

/  2  ,  2 
n/x  +  y 

begins  to  have  a  graph  whose  shape  is  more  like  that  of  the  function 

.  ,  2 
A  4*  X 

- — — .  Thus  numerical  results  assuming 

,.2.2 

4  +  x  -t  y 

/  2  2  ’ 
n/x  +  y  .  ■ 

may  be  approximated  by  a  polynomial  in  0  <  x  <  y-  are  not  very  good 

if  y  is  very  close  to  zero.  It  would  seem  that  a  similar  criticism' 

would  probably  apply  to  any  method  of  evaluating  the  singular  double 

integral  by  an  iterated  numerical  quadrature. 

The  other  possibility  is  to  set  up  a  single  two-dimensional 
formula  for  numerical  integration  of 

h  h 


i  =  55 


0  0  J2  .  ..2 


dxdy, 


(2.10) 


'x  +  y 

(where  F(x,y)  is  assumed  to  be  analytic  in  0  <  x  <  h,  0  <  y  <  h, 
and  where  F(0,0)  /  0) .  An  extremely  simple  symmetrical  formula  of 
this  type  is  of  the  form 


I  =  a  F(0,0)  +  b(F(0,h)  +  F(h,0))  +  c  F(h,h)  , 


(2.11) 


where  the  coefficients  a,b,  and  c  are  determined  (by  the  method  of 

undetermined  coefficients)  so  that  the  formula  is  exact  if  F  is 

any  polynomial  which  is  linear  in  both  the  x-andTy  directions.  That 

is,  formula  (2.11)  is  to  be  exact  if  F(x,y)  =  c.  +  c„x  +  c_y  +  c,xy. 

1  .  <  )  b 
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Since  the  last  funotion  is  linear,  it  is  only  necessary  to  require  that 
formula  (2.11)  be  exact  if  F(x,y)  =  1,  if  F(x,y)  =  x,  if  F(x,y)  =  y, 
and  if  F(x,y)  =  xy.  If  F(x,y)  =  1,  the  requirement  is 


+  2b  +  c  =  5  5  dxdy  =  2h  lnC/2  +  1) . 


0  0  /2  ,  2 
s/x  +  y 


(2.11a) 


If  F(x,y)  =  x,  one  obtains 

hb  +  he  =  J  5  *  y—  =  \  h2f/2-  1  +  In  (s/2  +  l)  ] . 

0  0  JT77  :  . 


(2.11b) 


If  F(x,y)  =  xy,  .  one  obtains 

h2c  =  55  -*?  -dxdy  =  |  h3(s/2  -  1).  (2.11c) 

0  0  /  2  .  2  J  .  • 
s/x  +  y 

If  F(x,y)  =  y,  because  of  symmetry  the  same  equation  is  obtained  as 
was  obtained  when  F(x,y)  =  x.  Solving  for  c,b,  eind  a, 

c  =  |  h(s/2  .-1) 

a  =  2b  =  j  h[3  ln(  1  +  s/2)  -s/2+1]. 

Hence,  formula  (2.1l)  may  be  conveniently  written 

I  =  h[. 276142375  F(h,h)  +  .371651200( F(0,h)  +  2F(0,0)  +  F(h,o))].  (2.12). 

In  a  similar  manner,  a  symmetrical  nine-point  formula  of  the  type 
I  =  a  F(0,0)  +  b[F(0,|)  +  F(|,0)  ]  +  c[F(0,h)  +  F(h,0)  ] 

+  d[F(|,h)  +  F(h,|)]  +  e  F(|,|)  +  g  F(h,h) 


(2.13) 
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may  be  obtained,  which  will  be  exact  for  any  polynomial  of  the  form 

2  2  2  2  2  2 
F(x,y)  =  cQ  +  c-jX  +  c2y  +  c^xy  +  c^x  +  c$y  +  c6x  y  +  c?xy  +  cgx  y  . 

The  coefficients  turn  out  to  be 

a  ~  jjj[ll  ln(s/2  +  1)  -s/2]  -  .276029863h 
.  _h_' 

30[13  In  (s/2  +1)  -  23  s/2  +  30]  =  .297698l57h 
c  =  ^[ln(s/2  +  1)  -  11  s/2  +  15]  -  .010834U7h 
d  =  j^[3  In  (s/5  +  1)  +  7  s/2  -  10]  =  .084787190h 
e  =  j^[24  In  (s/2  +  1)  +  56  s/2  T  80]  =  .678297519h 
’  g  =  9.  In  (s/2  +  1)  -s/2  +  10]  =  .021780805h. 

For  different  point  configurations  it  is  often  possible  to  derive 
formulas  of  this  type..  For  example,  the  seven -point  formula 

I  =  AF(0,0)  +  B[F(0,|)  +  F(|,0)]  +  C[F(0,h)  .+  F(h,0)]  +  DP(|,|)  +  EF(h,h) 
where 

•  A  '=  ^[5  ln(  1  +J2)  +  s/2  -  2]  =  .3l8423458h  - 

B  =  j^[4  ln(  1-  +s/2)  -  12  s/2  +  16]  =  .212910967h 

C  =  J^[ln(l  +  s/2)  -  3  s/2  +  4],  -  .053227742h 

D  =  ~[3  ln(l  +s/2)  +  7  s/2  -  10]  ■=  .847871899h 

■  E  =  3  ln(l  +  s/2)  +s/2  +  2] -i  . 064174400h , 

is  exact  if  F(x,y)  is  any  polynomial  of  the  form 

O  o  0  0 

F(x,y)  =  cQ  +  c1x  +  Cgy  +  c^xy  +  c^x  +  c^y  +  c^xy  +  c?x  y. 


(2.14) 
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From  the  way  the  coefficients  in  formulas  (2.13)  and  (2.14)  were 
obtained,  it  is  obvious  that  both  formulas  would  give  exact  results 
for  the  test  problem  of  Example  4.  A  new  example  is  therefore  considered t 

i  ±  e7x2+y 2 

Example  5 :  I,  =  5  5  . -  ■— dxdy.  (2.15) 

0  0  J? ~ 

The  integration  here  can  probably  not  be  performed  analytically. 

However,  upon  changing  to  polar  coordinates,  it  is  found  that 

I,  =  2  J  //4e*sec  9d0  -  £,  (2.16) 

0 

and  this  non-singular  integral  may  be  evaluated  quite  accurately  by 
ordinary  Gauss  formulas,  Simpson's  rule,  or  similar  formulas.  In  Table  2 
(on  the  next  page),  the  results  of  applying  various  types  of  formulas 
for  the  numerical  evaluation  of  the  integral  in  (2.15)  are  listed. 

Comments  on  the  results:  (l)  The  last  case  (in  which  the  singularity 

was  ignored)  was  included  merely  for  comparison.  When  ordinary  Gauss 

formulas  for  integration  of  non-singular  integrands  are  applied  like 

this  to  the  integration  of  non-negative  singular  integrands,  it  often 

happens  that  the  expected  poor  "answers"  are  much  too  small.  (2)  Since 
/  2+  2 

the  function  e^  y  is  strongly  concave  upward,  it  is  to  be  expected 
that  formula  (2.12)  (which  uses  only  the  four  comer  points  and  which 
assumes  the  function  is  linear  between  them) ,  will  give  results  which 
are  much  too  large.  The  results  using  this  four-point  formula  are 
worse  than  those,  obtained  by  using  the  sixteen-point  formula  which 
ignores  the  singularity.  This  illustrates  the  obvious  fact  that  a  low 
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5 -place  Value 

Error 

Correct  Value:  (obtained  by  accurate  evalua¬ 

tion  of  the  integral  in 
(2.16)) 

.509626 

4-point  process:  Cubature  formula  (2.12) 

.52274 

.01311 

7-point  process:  Cubature  formula  (2.14) 

.51021 

.00058 

9-point  process:  Iterated  3-point  Gaussian 
quadrature  assuming  logarithmic  singularity 
in  the  integration  with  respect  to  y;  (a 
set-up  like  formulas  (2.4)  and  (2.5)  with 
lower  limit  zero) 

.50865 

-  .00098 

9-point  process:  Cubature  formula  (2.13) 

.51061 

. 00118 

16-point  process:  Iterated  /-point  Gaussian 
quadrature  assuming  logarithmic  singularity 
in  the  integration  with  respect  to  y;  (a 
set-up  like  formulas  (2.4)  and  (2.5)  with, 
lower  limit  zero) 

.50940 

-  .00023 

16-point  process:  Cubature  formula  (II-3) 
in  Appendix  II 

.51011 

.00048 

16-point  process:  Iterated  4-point  ordinary 
Gaussian  quadrature  ignoring  the  presence  of 
the  singularity  completely 

.50373 

-  .00590 

Table  2  -  Results  of  Numerical  Integration  of  1^ 

order  formula  should  not  be  used  over  a  large  region.  Generally, 
formula  (2.12)  would  only  be  used  in  a  small  region  around  the  singu¬ 
larity;  ordinary  formulas  for  evaluation  of  non-singular  integrals  may 
be  applied  over  the  remainder  of  the  region.  These  remarks  apply  also 
to  somewhat  higher  order  formulas.  The  meaning  of  "small  region  about 
the  singularity"  will  depend  on  how  widely  F(x,y)  varies  in  the 
region.  Intuitively,  one  asks  the  question:  "May-  F(x,y)  be  approxi- 
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mated  reasonably  well  in  this  region  by  a  polynomial  of  the  degree 
for  which  the  formula  is  exact?"  (3)  The  fact  that  the  7-point 
cubature  formula  gave  better  results  than  the  9-point  cubature  formula 
probably  occurred  only  by  chance.  (4)  In  comparing  the  iterated 
quadrature  methods  used  in  this  example  with  the  cubature  methods 
of  formulas  (2.13)  and  (II-3)  of  Appendix  II,  it.  should  be  noted  that 
the  following  considerations  will  have  a  lot  to  do  with  the  accuracy 
of  the  results: 

(a)  As  mentioned  before,  the  Gaussian  iterated  quadrature 
method,  suffers  from  the  fact  that  near  one  end  of  the 
y -interval,  the  x  integration  may  be  considered  to  be 
"near  singular",  although  .we  go  ahead  and  use  ordinary 
Gaussian  quadrature  anyway. 

'  (b)  Formulas  (2.13)  and  (II-3)  are  "equally -spaced"  formulas, 

and  such  formulas  do  not  have  the  accuracy  that  generalized 
Gauss  formulas  possess}  that  is,  they  are  not  exact  for 
polynomials  of  as  high  a  degree  as  is  a  corresponding  Gauss 
type  formula  using  the  same  number  of  points. 

.For  this  particular-  example,  it  is  seen  that  the  iterated  quadrature 
methods  turned  out  to  be  slightly  better  than  the  cubature  method 
in  both  9-point  and  16-point  cases.  On  the  other  hand,  the  cubature 
formulas  (2.14)  and  (II-3)  are  somewhat  easier  to  use.  This  is 
especially  true  in  the  application  of  the  formulas  in  the  numerical 
solution  of  integral  equations;  in  this  case,  a  factor  in  the  inte¬ 
grand  contains  the  unknown  variable,  and  it  is  usually  much  more 
convenient  to  use  an  equally -spaced  formula  than  one  with  any  other 
spacing. 
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§3.  Singular  integrals  of  the  type 

1=5  5  In  Jx2  +  y2  F(x,y)dxdy  (3.1) 

0  0 

Here  it  is  assumed  that  F(x,y)  is  analytic  in  the  region 

<h  <  \J2 

and  that  F(0,0)  ^  0.  In  a  manner  analogous  to  that  used  in  paragraph 
2,  equally -spaced  cubature  formulas  may  be  derived.  In  Appendix  III, 
four-point,  nine-point,  and  sixteen-point  formulas  of  this  type  are 
listed.  It  is  to  be  noted  that  the  coefficients  in  these  formulas 
will  always  be  negative  if  h  is  in  the  range  [0,|-  s/2]  which  was 
specified,  (intuitively,  one  wishes  for  practical  use  to  rule  out 
any  formulas  which  have  some  negative  and  some  positive  coefficients.) 

A  good  example  of  a  problem  of  type  (3.1)  in  which  F(x,y)  is 
not  a  polynomial  and  for  which  the  exact  solution  is  known  has  eluded 
the  author.  In  the  absence  of  such  a  test  problem,  the  following  problem 
was  considered: 

h  h  /  p  o’  o 

Example  6:  J  =  5  5  In  Jx  +  y  cos(5xy  )dxdy  (3*2) 

'00 

For  various  values  of  h,  the  formulas  of  Appendix  III  were  tried 
on  this  problem.  Also,  high  order  iterated  Gaussian  quadrature  (com¬ 
pletely  ignoring  the  singularity)  was  tried.  The  results  of  evaluating 
-J  by  .  various  formulas  are  given  in  Table  3. 
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h  =  0.1 

h  =  0.2 

h  =  0.3 

h  =  0.5 

h  =  1.0 

4-point  formula 

(III-D 

.02670606 

.079065 

.1412 

.256 

.36 

9 -point  formula 
(II I -2) 

.02670611 

.07909542 

.141421 

.2635 

.3855 

16-point  formula 
(III -3) 

.02670611 

.07909538 

. 141419 

.263395 

•  3.900 

9 -point  formula 
(ignoring  singularity) 

.02673 

.0792 

i  1416 

.264 

•  394 

64-point  formula 
(ignoring  singularity) 

.026707 

.079099 

.14143 

.263455 

.390973 

256-point  formula 
(ignoring  singularity) 

.02670617 

.07909569 

.141421 

.263437 

.390900 

h  ^  y  2  2  2 

Table  3*  Results  in  numerical  integration  of  -  5  5  ln<Jx  +y  cos(5xy  )dxdy 

0  0 

Comments  on  the  results t  (1)  Although  the  exact  results  are  not  known, 
enough  work  has  been  done  with  this  problem  so  that  it  is  felt  that  we 
at  least  know  when  the  ’’answers”  are  obviously  quite  a  bit  off.  For 
example,  then,  when  it  was  felt  that  the  "answer"  was  certainly  good 
to  no  more  than  three  significant  figures,  only  this  many  figures,  were 
listed  in  the  table.  (2)  One  notes  that  the  three  cubature  formulas 
(111-1,2,3)  give  essentially  the  same  results  for  h  =  0.1.  For  larger 
values  of  h,  the  four-point  formula  soon  begins  to  give  noticeably 
poorer  results.  (The  function  F(x,y)  =  cos(5xy  )  cannot  be  approxi- 

*  X 

mated  very  well  in  an  interval  0  <  <  h  if  h  is  more  than  0.1  or 

y 

0.2.)  However,  the  nine-point  formula  and  the  sixteen-point  formula 
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give  quite  comparable  results  for  h  =  0.3*  For  h  =  1.0,  none  of 

the  three  cubature  formulas  should  have  been  used.  They  were  derived 

assuming  h  <  |-s/2,  and  when  one  uses  them  for  h  =  1.0,  the  coefficients 

in  the  formula  are  not  of  one  sign}  of  course,  the  weighting  function 
I~ 2  2 

In  sfx  +  y  is  not  always  non-positive  either  if  in  part  of  the  region 
x*  +  y  >1,  (which  will  happen  if  h  >  ^-s/2)  •  (3)  The  formulas  used 
to  obtain  the  results  in  the  last  three  rows  of  'the  Table  were  ordinary 
Gaussian  quadrature  formulas  applied  to  both  x  and  y  integrations. 

For  example,  for  the  last  line,  sixteen-point  Gauss  formulas  were  used 
for  each  quadrature.  Instead  of  our  problem,. if  one  were  trying  to 
evaluate 

h  h 

j  5  c os ( 5 xy  )dxdy ,  . 

0  0  ‘  . 

even  the  nine-point  formulas  would  give  excellent'  results  for  the 
smaller  values  of  h)  and  the  64-point  and • 256-point  formulas  would 
have  given  excellent  results-  for  all  values  of  h  listed.  But  when 
one  tries  to  apply  these  formulas  when  the  integrand  has  a  singularity, 
the  results  (as  expected)  are  not  good.  For  example,  for  h  =  0.3, 
the  256-point  formula  ignoring  the.  singularity  gives  the  same  results 
as  the  nine-point  cubature  'formula  (III-2).  ■ 

§4. ■  It  is  hoped  that  the  preceding  examples  will  cause  programmers 
to  look  carefully  at  any  integral  they  wish  to  evaluate  before  deciding 
what  integration  formula  to  use  on  the  computer. 
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Formulas  similar  to  those  listed  in  the  three  appendices  may  be 
derived  for  many  other  specific  problems.  It  is  certainly  possible 
to  set  up  a  machine  program  for  deriving  such  formulas.  Such  a  program 
would,  however,  usually  require  multiple  precision  arithmetic.  For 
example,  for  the  sixteen-point  cubature  formula  (II-3) ,  it  was  necessary 
to  solve  ten  linear  equations  for  the  ten  unknown  coefficients  in  the 
formula,  and  these  equations  are  quite  ill-conditioned.  For  the  four- 
point  formula  (il-l) ,  the  actual  equations  were  listed  earlier  in  this 
document  (equations  (2.11a)  , (2.11b) ,  and  (2.11c)).  This  particular 
system  is  already  in  triangular  form,  but  it  will  not  be  true  in  general 

From  equations  (2.11a)  -'(2.11c),  one  sees  that  another  difficulty 
in  deriving  formulas  for  specific  integration  problems  might  very  well 
be  in  the  evaluation  of  the  integrals  analogous  to  those  on  the  right 
in  equations  (2.11a)  r  (2.11c). 
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Appendix  I 

Table  of  Abscissas  and  Coefficients  in  the  Generalized  Gaussian  Formula 

h  n 

5  In  x  F(x)dx  =  Z  H,F(x.) 

0  i=l  1  1 


A  description  of  how  to  derive  these  formulas  is  given,  for 
example,  in  Kopal1  or  Mineur8.  The  function  F(x)  is  assumed  to 
have  2n  continuous  derivatives  in  [0,h],  where  h  <  1.  The  fact 
that  the  integral  only  exists  as  an  improper  integral  does  not  have 
any  effect  on  the  method  of  derivation',  but  it  does  limit  the  choice 
of  methods  of  deriving  simple  error  expressions. 

The  basic  idea  is  that  the  formula 
h  n 

i  In  x  F(x)dx  =  S  H.  F(x. )  +  R  (i) 

0  i=l  11  n 

is  derived  so  that  the  remainder  term  is  identically  zero  if 

F(x)  is  any  polynomial  of  degree  no  more  than  2n  -  1.  In  the  table 
below,  the  proper  values  of  the  abscissas  and  the  corresponding 

coefficients  are  given  for  various  fixed  values  of  h  and  n. 


The  Hermite  representation  for  F(x)  in  the  interval  [0,h] 
using  the  base  points  x-^,x,,,...,xn  is  (see  Steffensen3) 

F(x)  =  iZi[t1(x)]2fF(xi)[l  -  2(x  -  x1)^(xi)]  +(x  -  xi)F'(xi)} 


.  F(2n)  U(x)l  £r  ,  x2 

+  (2tfi  £ (x  -  Xj) 


(ii) 
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where 


^U) 


n 


=  IT 


k=l 

(tfi) 


(x  -  \\ 

<xi  - 


and  where  0  <  5(x)  <  h. 


Since  the  stun  in  equation  (ii)  (which  will  be  designated  as 
^2n  merely  a  polynomial  of  degree  2n  -  1,  and  since  F(x) 

is  continuous,  it  follows  that  F^^[£(x)]  is  a  continuous  function 
of  x.  Substitution  of  equation  (ii)  into  (i)  then  gives  an  expression 

for  V 

h  h  P(2n) r  rf  v\ l  n  O 

Rn  =  l  ln  x  *2n-l(x)dx  +  l  ln  X  Wl  JT  (x  "  Xj)  dX 

n  n  F^2n^[^(x.)]  n  2 

"i^i  HiQ2n-i(Xi^  ■  i21  Hi  ran  .  JJ^i  ’  xj  * 

The  product  in  the  last  term  is  zeroj  the  first  integral  is  equal 
to  the  first  sum  because  the  formula  was  originally  derived  to  be  exact 
for  polynomials  of  degree  (2n  -  l)  or  less.  Thus 


h 

=  5  ln  x 
0 


p(  2n)  [  £(x)  ] 

(2571 


JT  (x  -  x.)  dx. 
j=l  3 


Since  p(2n)[£(x)]  iB  continuous  and  since  the  other  factor  in  the 

integrand  is  always  negative,  one  may  use  the  mean  value  theorems 

applied  to  the  improper  integral  to  obtain  the  formula 

p(2n)m  h  n  ? 

Rn  =  (2n)l  l  ln  x  TT  U  "  XJ)  dX 
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vhere  0  £  C  £  h.  The  integration  here  may  be  carried  out  in  closed 
form,  and  for  any  particular  values  of  h  and  n  a  numerical  value 
of  the  coefficient  of  p(2n)(£)  could  be  obtained  if  desired. 


Abscissas  and  coefficients  for 


the  formulas  are  given  in  the 


following  tables. 


h 

5  In  x  F(x)dx 
0 


Two-point  formulas 


xi 

Hi 

.01.7  065  3648 

CO 
CO 
•1— 1 

• 

1 

936 

h'  = 

.1 

.076  340  9457 

-.141 

321 

.032  526  740 

-309 

155 

h  = 

.2 

.151  048  729 

-.212 

732 

.039  827  506 

-.358 

810 

h  = 

.25 

.187  818  171 

-.237 

762 

.046  866  292 

-.403 

476 

h  = 

•  3 

.224  159  419 

-.257 

715 

.060  185  480 

-.480 

899 

h  = 

.4 

.295  361  084 

-.285 

616 

.072  495  937- 

-.545 

589 

h  = 

.5 

.364  106  635 

-.300 

984 

.098  332  529 

-.664 

069 

h  = 

.75 

.517  865  507 

-.301 

692 

n 

=  2  H  F(x.) 

i=l  ■  1 


Three-point  formulas 


xi 

Hi 

512 

.009 

038 

5371 

-.117 

704 

182 

998 

.046 

848 

5242 

-.139 

364 

134 

.087 

717' 

3234 

-.073 

190 

193 

272 

.017 

311 

2632 

-.197 

105 

695 

310 

.092 

032 

2473 

■  -.217 

367 

372 

.174 

760 

3835 

-.107 

414 

515 

887 

.021 

252 

272 

-.230 

944 

257 

703 

.114 

098 

280 

-.247 

084 

847 

.218 

033 

359 

-.118 

544 

486 

458 

.025 

076 

930  ‘ 

-.261 

970 

538 

384 

.135 

810 

533 

-.272 

407 

593 

.261 

119 

549 

-.126 

813' 

710 

476 

.032 

395 

246 

—317 

275 

288 

817 

.178 

113 

942 

-.312 

698 

376 

.346 

625 

133 

-.136 

542 

629 

256 

.039 

281 

165  ■ 

-.365 

307 

501 

334 

.218 

757 

540 

-.342 

191 

529 

•  430 

968 

225 

-.139 

074 

559 

402 

.054 

430 

334 

-.460 

263 

416 

152 

.310 

616 

188 

-.382 

620 

062 

.631 

807 

628 

-.122 

878 

077 

.112  006  806 
h  =  1.0  .602  276  908 


-.718  539  319 
—  281  460  681 


063  890  793  -.513  404  552 
368  997  064  -.391  980  041 
766  880  304  —  094  615  407 
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h 

S  In  x  F(x)dx 
0 


Four-point 

Formulas 

Xi 

Hi 

.005  579  3787 

-.080  432 

.030  456  3106 

-.112  913 

h  = 

.1 

.065  083  2445 

-.092  721 

.092  564  6649 

-.044  191 

.010  730  0607 

-.136  781 

.059  732  9846 

-.180  961 

h  = 

.2 

.129  046  5053 

-.140  154 

.184  800  1041 

-.063  989 

.013  199  0687 

-.161  254 

.074  022  3111 

-.208  196 

h  = 

.25 

.160  649  9113 

-.156  984 

.230  795  6060 

-.070  137 

.015  605  5281 

-.183  944 

.088  086  1632 

-.232  212 

h  = 

.3 

.191  986  7991 

-.170  544 

.276  698  6323 

-.074  489 

.020  242  153 

-.225  038 

.115  524  836 

-.272  635 

h  = 

.4 

.253  774  725 

-.189  937 

.368  171  418 

-.078  904 

.024  650  669 

-.261  493 

.141  985  961 

-.305  065 

h  = 

•  5 

.314  166  312 

-.201  254 

;459  053  527 

-.078  760 

.034  597  958 

-.336  501 

.202  794  754 

-.360  747 

h  = 

.75 

.455  812  835 

-.204  766 

.680  823  438 

-.063  747 

.041  448  480 

-.383  464 

.245  274  914 

-.386  875 

h  = 

1.0 

.556  165  454 

-.190  435 

.848  98  2  395 

-.039  225 

n 

Z  H,F(x.) 
i=l  1  1 


Five -point  Formulas 


xi 

Hi 

5170 

.003 

785 

3873 

-.058 

679 

2293 

0245 

.021 

146 

0155 

-.090 

010 

3886 

2481 

.048 

016 

9933 

-.087 

763 

4421 

7197 

.075 

741 

5540 

-.064 

349 

1963 

.095 

035 

7147 

-.029 

456 

2529 

961 

.007 

303 

6269 

-.100 

898 

343 

542 

.041 

473 

7959 

-.147 

093 

364 

911 

.094 

989 

2687 

-.136 

215 

523 

169 

.150 

759 

6325 

-.095 

367 

634 

.189 

888 

8128 

-.042 

312 

718 

261 

.008 

997 

7473 

-.119 

469 

913 

736 

.051 

405 

8244 

-.170 

641 

071 

793 

.118 

143 

6743 

-.154 

482 

727 

800 

.188 

016 

1836 

-.105 

790 

569 

.237 

247 

4844 

-.046 

189 

311 

963 

.010 

653 

955 

-.136 

818 

022 

713 

.061 

191 

357 

-.191 

814 

948 

383 

.141 

073 

659 

-.169 

924 

693 

782 

.225 

090 

424 

-.113 

791 

218 

.284 

554 

707 

-.048 

842 

960 

822  ■ 

.013 

860 

195 

-.168 

565 

160 

015 

.080 

323 

520 

-.228 

532 

902 

617 

.186 

217 

509 

-.194 

139 

340 

839 

.298 

616 

333 

-.124 

050 

550 

.378 

983 

205 

-.051 

228 

341 

609 

.016 

929 

963 

-.197 

115 

290 

815 

.098 

846 

742 

-.259 

304 

126 

116 

.230 

289 

847 

-.211 

378 

918 

050 

.371 

122 

736 

-.128 

274 

873 

.473 

060 

262 

-.050 

500 

383 

008 

.023 

968 

326 

-.257 

348 

518 

058 

.141 

941 

301 

-.316 

979 

130 

138 

.334 

012 

731 

-.232 

815 

086 

351 

.544 

908 

302 

-.120 

066 

808 

.705 

185 

931 

-.038 

552 

012 

068 

.029 

134 

472 

-.297 

893 

472 

318 

.173 

977 

213 

-.349 

776 

227 

127 

.431 

702 

520 

-.234 

488 

290 

487 

.677 

314 

175 

-.098 

930 

459 

.894 

771 

361 

-.018 

911 

552 
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I 


h  n 

5  In  x  F(x)dx  =  Z  H.  F(x. ) 
0  i=l 


h  =  .1 


h  =  .2 


h  =  .25' 


h  =  .3 


h  =  .4 


h  =  .5 


Slx-polnt  Formulas 


xi 

Hi 

.002  737 

2522 

-.044 

867 

1463 

.015  469 

6555 

—  072 

651 

9894 

.036  303 

1065 

-.077 

680 

9568 

.060  487 

2300 

-.067 

222 

4997 

.08  2  296 

5498 

-.046 

834 

4688 

.096  456 

8465 

-.021 

001 

4482 

.005  295 

0647 

-.077 

799 

6177 

.030  357 

5278 

-.120 

458 

4218 

.071  748 

1227 

-.123 

092 

4025 

.120  160 

3044 

-.102 

007 

6953 

.164  114 

1637 

-.068 

577 

9746 

.192  802 

6065 

-.030 

Oil 

4705 

.006  530 

9449 

-.092 

422 

0578 

.037  6a 

3068 

-.140 

584 

7603 

.089  209 

6762 

-.140 

908 

2984 

.149  727 

1075 

-.114 

469 

7378 

.204  852 

8015 

-.075 

514 

5934 

.240  934 

2157 

-.032 

674 

1424 

.007  7U 

8903 

-.106 

154 

737 

.044  825 

428 

-.158 

911 

943 

.106  499 

312 

-.156 

404 

554 

.179  106 

115 

-.124 

589 

891 

.245  466 

855 

-.080 

676 

771 

.289  034 

540 

-.034 

453 

945 

.010  094 

293 

-.131 

473 

543 

.058  898 

892 

-.191 

288 

667 

.140  543 

490 

-.181 

914 

363 

.237  .247 

089 

-.139 

293 

045 

.326  262 

632 

-.086 

642 

714 

.385  122 

046 

-.035 

903 

961 

.012  357 

848 

-154 

463 

257 

.072  569 

130 

-.219 

128 

433 

.173  810 

761 

-.201 

672 

708 

.294  428 

876 

-.148 

195 

050 

.406  333 

465 

-.088 

009 

317 

.481  007 

530 

-.035 

104 

825 

Seven-point  Formulas 


xi 

Hi 

.002 

072 

1068 

-.035 

527 

6951 

.011 

784 

1174 

-.059 

663 

9342 

.028 

206 

4839 

-.067 

568 

1217 

.048 

553 

9795 

-.063 

982 

2229 

.069 

244 

2584 

-.032 

329 

9697 

.086 

559 

5443 

-.035 

470 

0086 

.097 

346 

58  00 

-.015 

716 

5571 

.004 

016 

8810 

-.062 

016 

4469 

.023 

ia 

2079 

-.100 

038 

5777 

.055 

729 

0893 

-.108 

818 

0112 

.096 

347 

1971 

-.099 

062 

9590 

.137 

872 

9565 

-.078 

139 

9389 

.172 

790 

5153 

-.051 

433 

2069 

.194 

620 

7986 

-.022 

378 

4418 

.004 

959 

0934 

-.073 

862 

5772 

.028 

704 

7253 

-.117 

287 

7392 

.069 

288 

9353 

-.125 

452 

0207 

.120 

001 

7469 

-.112 

217 

0601 

.171 

977 

8038 

-.087 

005 

2591 

.215 

788 

0211 

-.056 

429 

8697 

.243 

231 

0509 

-.024 

319 

0640 

.005 

883 

9101 

-.065 

033 

5373 

.034 

197 

4598 

-133 

134 

9773 

.082 

718 

5014 

-.140 

190 

7592 

.143 

492 

1782 

-.123 

289 

2070 

.205 

933 

9771 

-.093 

948 

6461 

.258 

698 

4269 

-.060 

001 

0805 

.291 

820 

9534 

-.025 

593 

6338 

.007 

685 

2608 

-.105 

745 

899 

.044 

975 

476 

-.161 

490 

108 

.109 

178 

862 

-.165 

181 

407 

.189 

947 

832 

-.140 

512 

231 

.273 

347 

751 

-.103 

282 

880 

.344 

213 

964 

-.063 

752 

179 

oo 

$ 

927 

180 

-.026 

551 

589 

.009 

425 

1483 

-.124 

689 

545 

•  055 

472 

574 

-.186 

293 

354 

.135 

071 

476 

-.185 

450 

850 

.235 

614 

042 

-.152 

577 

010 

.339 

964 

610 

-.107 

834 

450 

.429 

209 

778 

—  063 

915 

195 

.485 

902 

001 

-.025 

813 

186 

L 
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h  n 

J  In  x  F(x)dx  =  Z  H.  F(x. ) 
0  i=l  1  1 


Six-point  Formulas  Seven-point  Formulas 


xi 

Hi 

Xi 

Hi 

.017 

606 

007 

-.203 

816 

913 

.013 

492 

725 

—  165 

885 

575 

.104 

662 

873 

-.273 

867 

828 

.080 

287 

312 

—236 

580 

830 

.252 

541 

235 

-.232 

953 

197 

.196 

665 

459 

-.221 

087 

185 

.431 

126 

864 

-.152 

443 

361 

•  344 

963 

158 

-166 

619 

048 

.600 

829 

037 

-.076 

930 

220 

•  500 

964 

878 

-.104 

280 

334 

.718 

785 

844 

-.025 

750 

035 

.637 

433 

036 

-.052 

901 

605 

.727 

073 

965 

-.018 

406 

978 

.021 

634 

006 

-.238 

763 

663 

.016 

719 

355 

-.196 

169 

389 

.129 

583 

391 

-.308 

286 

573 

.100 

185 

678 

-.270 

302 

644 

.314 

020 

450 

-.245 

317 

427 

.246 

294 

246 

-239 

681 

873 

.538 

657 

217 

-.142 

008 

757 

.433 

463 

493 

—  165 

775 

775 

.756 

915 

337 

-.055 

454 

622 

.632 

350 

988 

-.088 

943 

227 

.922 

668 

851 

-.010 

168 

959 

.811 

118 

627 

-033 

194 

304 

.940 

848 

167 

-.005 

932 

7870 

Eight-point  Formulas 

.001 

623 

6190 

n 

o 

1 

901 

5876 

.009 

265 

6696 

-.049 

825 

3739 

.022 

462 

9083 

-.058 

678 

5672 

.039 

487 

9968 

-.058 

818 

6615 

.058 

025 

2058 

-.052 

495 

3423 

.075 

504 

4191 

-.041 

603 

7363 

.089 

469 

1394 

-.027 

737 

3689 

.097 

939 

6237 

-.012 

197 

8715 

.003 

153 

0217 

-.050 

725 

2236 

.018 

208 

5251 

-.084 

294 

1729 

•  044 

381 

3184 

-.095 

731 

6168 

.078 

310 

8566 

—  092 

590 

0626 

.115 

415 

5869 

-.079 

862 

3292 

.150 

543 

2585 

-.061 

395 

6310 

.178 

704 

9274 

-.039 

965 

9145 

.195 

829 

5772 

-.017 

322 

6319 

.003 

895 

6151 

-.060 

541 

1921 

.022 

594 

4803 

-099 

185 

5690 

.055 

184 

0243 

-.110 

975 

3784 

.097 

517 

2457 

-.105 

675 

7286 

.143 

902 

4230 

—  089 

725 

6908 

.187 

901 

3513 

-.067 

964 

3671 

.223 

238 

4836 

-.043 

706 

0875 

.244 

756 

1390 

-.018 

799 

5766 

29 


h  n 

5  In  x  F(x)dx  =  2  H,F(x.) 

0  i=l 


Eight-point  Formulas 


Xi 

Hi 

.004 

625 

5195 

-.069 

827 

8231 

.026 

928 

2927 

-.112 

957 

6094 

.065 

886 

9117 

-.124 

658 

4226 

.116 

587 

3281 

-.116 

957 

3535 

.172 

245 

0299 

-.097 

773 

8166 

.225 

143 

4404 

-.072 

948 

1768 

.267 

7C9 

4464 

-.046 

311 

9302 

.293 

668 

7565 

-.019 

756 

7089 

.006 

050 

2310 

-.087 

123 

271 

.035 

444 

175 

-.137 

831 

802 

.066 

990 

833 

-.148 

324 

262 

.154 

299 

396 

-.135 

253 

690 

.228 

458 

891 

-.109 

574 

963 

.299 

232 

903 

-.079 

155 

968 

.356 

430 

573 

-.048 

822 

802 

.391 

443 

600 

-.020 

429 

534 

•  007 

430 

4341 

-.103 

031 

774 

.043 

756 

172 

-.159 

860 

455 

.107 

670 

933 

-.168 

086 

490 

.191 

382 

639 

-.149 

070 

658 

.283 

944 

674 

-.116 

869 

081 

.372 

678 

364 

-.081 

404 

215 

.444 

772 

010 

-.048 

471 

625 

.489 

128 

485 

-.019 

779 

292 

.010 

677 

906 

-.137 

975 

753 

.063 

512 

116 

-.205 

495 

866 

.157 

081 

786 

-.204 

954 

207 

.280 

413 

715 

-.169 

505 

437 

.417 

959 

377 

-.121 

017 

820 

.551 

576 

813 

-.074 

576 

282 

.662 

404 

789 

-.038 

421 

023 

.732 

473 

320 

-.013 

815 

167 

.013 

320 

244 

-.164 

416 

605 

.079 

750 

429 

-237 

525 

610 

.197 

871 

029 

-.226 

841 

984 

.354 

153 

994 

-.175 

754 

079 

.529 

458 

575 

-.112 

924 

030 

.701 

814 

530 

-.057 

872 

211 

.849 

379 

320 

-.020 

979 

074 

.953 

326 

450 

-.003 

686 

4071 
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Appendix  II 


h  h 

Cubature  Formulas  for  1=5  5 

0  0 


F(*,y) 


[~2~,  2 
+  y 


dxdy  where 


0  <  h  <  ~J2 


Four  Point  Formula 

I  =  a],F(0,0)  +  a2[F(h,0)  +  F(0,h)]  +  a3F(h,h)  where 

a1  =  j[3  ln(l  +J2)  -  Jz  +  1]  =  .743  302  400h 

a2  =  |[3  ln(l  +J2)  -  Jz  +  1]  =  .371  651  200h  (H-l) 

a3  =  y<  J2  -  1)  =  .276  142  375h 


The  above  formula  is  exact  if  ?(x,y)  is  any  polynomial  of  the  form 
c0  +  clx  +  ^  +  c3Xy* 


Nine  Point  Formula 

I  =  *.^(0,0)  +  a2[F(|,0)  +  F(0,|)]  +  a3[F(h,0)  +  F(0,h)] 

+  a4F(|,|)  +  a5[F(h,|)  +  F(|,h)]  +  a6F(h,h) 

where  (when  L  =  ln(l  +  «/2)  and  when  S  =  and  when  K  =  jr) 

a3  =  11L  -  S  =  .2760  298 63h 

a2  =  13L  -  23S  +  6K  =  .2976  98l57h 

a3  =  L  -  US  +  3K  =  .0106  34147h  (II-2) 

a4  =  24L  +  56S  -  16k  =  .6782  97519h 

a5  =  3L  +  7S  -  2K  =  .0847  87190h 

a6  =  -  9L  -  S  +  2K  =  .0217  80805h 

The  above  formula  is  exact  if  F(x,y)  is  any  polynomial  of  the  form 
c0  +  Clx  +  c2x2  +  c3y  +  c4xy  +  c^y  +  c6y2  +  c?xy2  +  CgX2y2. 
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Appendix  III 
b  h  i  2  o 

Cubature  Formulas  for  1  =  5  J  F(x,y)  ln^x  +  y  dxdy  where 

0  0 


0  <  h  <  \ 


Four  Point  Formula 

I  =  a^FCOjO)  +  a2[F(h,0)  +  F(0,h)]  +  a^F(h,h)  where 
2 

al  =  ln  h  +  4  In  2  +  4tt  -  25]  =  h2[|  In  h  -  .201  271  680] 

2 

a2  =  !M12  ln  h  +  4  in  2  +  4n  -  19]  =  h2[j  ln  h  -  .076  271  680]  (III-1) 

2 

a3  =  W^-12  ln  h  +  12  ln  2  -  9]  =  h2[|  ln  h  -  .014  213  205] 

The  above  formula  is  exact  if  F(x,y)  is  any  polynomial  of  the  form 
C0  +  °1X  +  C2y  +  C3xy< 

Nine  Point  Formula 

I  =  a1F(0,0)  +  a2[F(|,0)  +  F(0,|)]  +  a.,[F(h,0)  +  F(0,h)] 

+  a4F(|,|)  +  a5[F(h,|)  +  F(|,h)]  +  a^h.h) 

where  (when  H*  =  h  h  and  L  =  h  -j^-2  and  P  =  and  K  = 

&1  =  H*  +  7L  +  2P  -  79K  =  h2[^-  ln  h  -  .065  313  206] 

a2  =  4H*  -  20L  +  20P  -  140K  =  h2[g-  ln  h  -  .096  927  873] 

a3  =  H*  -  17L  +  8P  -  K  =  h2[jg-  ln  h  +  .002  960  3808]  (III-2) 

a4  =  16H*  +  112L  +  32P  -  624K  =  h2[|ln  h  -  .156  122  407] 

a$  =  4H*  +  28L  +  8P  -  116K  =  h2[i-  ln  h  +  .016  524  954] 

a6  =  H*  -  11L  -  16P  +  137K  =  h2[^r  ln  h  +  .008  292  4433] 

The  above  formula  is  exact  is  F(x,y)  Is  any  polynomial  of  the  form 

c0  +  clx  +  c2x2  +  c3y  +  c4Xy  +  c5x2y  +  c6y2  +  c7Xy2  +  c8x2y2* 


Sixteen  Point  Formula 

I  =  a^O.O)  +  a2[F(j,0)  +  F(0,j)]  +  a3[F(y,0)  +  F(0,^)] 

+  a4[F(h,0)  +  F(0,h)]  +  a5F(j,|)  +  «6[Ky»j)  + 

+  a?[F(h,  j)  +  F(j,h)]  +  agF(y,y)  +  ag[F(h,^-)  +  F(y,h)] 

+  a1(/(h*h) 

where  (when  L  -  h  JS  „d  S  =  and  K  = 

a^  =  28L  -  292S  +  236K  =  .171  814  675h 

a2  =  75L  -  1413S  +  657K  =  .209  487  987h 

a3  =  -  114L  -  1170S  +  1926K  =  .026  000  525h  (II-3) 

a4  =  -  49L  -  65S  +  541K  =  .024  744  850h 

a5  =  5184S  -  2916K  =  .446  215  211h 

a6  =  405L  -  3483S  -  1053K  =  .138  207  298h 

a?  =  18 OL  -  1548S  -  468k  =  .061  425  466h 

ag  =  -  648L  +  9720S  -  648K  =  .135  840  492h 

ag  =  -  63L  +  1233S  -  225K  =  .037  996  448h 

a10  =  72L  "  1720S  +  572K  =  ,013  151  648h 

The  above  formula  is  exact  if  F(x,y)  is  any  polynomial  of  the  form 
c0  +  clx  +  c2x2  +  c3x3  +  c4y  +  c5xy  +  c6x2y  +  c7x3y  +  c8y2  +  c9xy2 

+  ciox2y2  +  cnx3y2  +  ci2y3  +  ci3xy3  +  cux2y3  +  cis^y3- 
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Appendix  III  (continued) 


Cubature  Formulas  for 


h  h 

I  =  $  S  F(x,y) 
0  0 


In 


J? 


+  y  dxdy 


where 


0  <  h  <  |  s/2 


Sixteen  Point  Formulas 

I  =  a1F(0,0)  +  a2[F(j,0)  +  F(0,j)]  +  a3[F(yfO)  +  F(0,y)] 

+  a4[F(h,0)  +  F(0,h))  +  a$F(|,j)  +  *6[F(y,j)  +  F(j,y)] 

+  a?[F(h,  j)  +  F(j,h)]  +  agF(y,y)  +  ag[F(h,^)  +  F(y,h)]  +  a 

2  2  2 
where  (when  H*  =  — and  L  =  and  p  =  and  K  = 

a^  =  H*  -  13L  -  4P  -  967K  =  hZ[g£  In  h  -  .030  830  973] 

a2  =  3H*  +  27L  +  81P  -  12207K  =  h* I 2[^-  In  h  -  .064  455  221] 

a^  =  3H*  -  351L  -  108P  +  18411K  =  h2[g£  In  h  -  .004  311  2601] 

a4  =  H*  -  139L  -  67P  +  9659K  =  h2[^-  In  h  -  .002  970  4856] 

a5  =  9H*  +  270L  -  135P  -  1215K  =  h2[^  In  h  -  .163  651  202] 

a6  =  9H*  +  513L  +  59 4P  -  72333K  =  h2[^  In  h  -  .023  045  066] 

a?  =  3H*  +  198L  +  279P  -  328  53K  =  h2[^-  In  h  -  .007  683  7845] 

ag  =  9H*  +  27L  -  864P  +  86913K  =  h2[£j  In  h  +  .012  148  911] 

ag  =  3H*  +  63L  -  126P  +  11697K  =  h2[^  In  h  +  .007  951  6953] 

a10  =  H*  "  +  117P  -  10119K  =  h2[^-  In  h  +  .003  333  2599] 


1QF(h,h) 

h2 

h  \ 
537So' 


(III -3) 


The  above  formula  is  exact  if  F(x,y)  is  any  polynomial  of  the  form 

2  3  2  3  2  2 

c0  +  clx  +  C2X  +  C3X  +  °4y  +  c5xy  +  C6X  y  +  C7X  y  +  c8y  +  c9xy 
22  32  3  3  23  33 

+  c10x  y  +  cllx'^y  +  c12y  +  C13xy  +  c14x  y  +  c15^y  * 


34 
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